Neuroinflammation, Energy and Sphingolipid Metabolism Biomarkers Are Revealed by Metabolic Modeling of Autistic Brains

Autism spectrum disorders (ASD) are a heterogeneous group of neurodevelopmental disorders generally characterized by repetitive behaviors and difficulties in communication and social behavior. Despite its heterogeneous nature, several metabolic dysregulations are prevalent in individuals with ASD. This work aims to understand ASD brain metabolism by constructing an ASD-specific prefrontal cortex genome-scale metabolic model (GEM) using transcriptomics data to decipher novel neuroinflammatory biomarkers. The healthy and ASD-specific models are compared via uniform sampling to identify ASD-exclusive metabolic features. Noticeably, the results of our simulations and those found in the literature are comparable, supporting the accuracy of our reconstructed ASD model. We identified that several oxidative stress, mitochondrial dysfunction, and inflammatory markers are elevated in ASD. While oxidative phosphorylation fluxes were similar for healthy and ASD-specific models, and the fluxes through the pathway were nearly undisturbed, the tricarboxylic acid (TCA) fluxes indicated disruptions in the pathway. Similarly, the secretions of mitochondrial dysfunction markers such as pyruvate are found to be higher, as well as the activities of oxidative stress marker enzymes like alanine and aspartate aminotransferases (ALT and AST) and glutathione-disulfide reductase (GSR). We also detected abnormalities in the sphingolipid metabolism, which has been implicated in many inflammatory and immune processes, but its relationship with ASD has not been thoroughly explored in the existing literature. We suggest that important sphingolipid metabolites, such as sphingosine-1-phosphate (S1P), ceramide, and glucosylceramide, may be promising biomarkers for the diagnosis of ASD and provide an opportunity for the adoption of early intervention for young children.


Introduction
Autism spectrum disorder (ASD) is a neurodevelopmental disorder that shows up in early childhood and is a lifelong condition. Although the heterogeneous nature of this disorder creates difficulties when diagnosing people, most ASD patients lack verbal and nonverbal communication skills, insist on the same routine, and show restricted and repetitive behaviors [1]. The prevalence of ASD among children varies worldwide, ranging from 0.01% to 4.36%, with higher prevalence values in high-income countries. Methodological, cultural, and environmental variables, as well as access to mental healthcare and research funding, may all play a role in variations in prevalence values. Nevertheless, two conclusions were prominently similar in all reports: the prevalence of ASD increases year by year, and males are, on average, four times more susceptible than females [2].
A number of studies emphasized that various brain parts, such as the frontal lobes, amygdala, and cerebellum, may show structural and functional differences in ASD and can take part in the pathology of the disorder [3]. For instance, the cortical thickness of the frontal cortex was found to be increased in the frontal cortex but decreased in the

ASD-Specific Prefrontal Cortex GEM Construction
The postmortem prefrontal cortex transcriptome data for autistic people and neurotypical controls aged between 2 and 56 were acquired from the Gene Expression Omnibus database with the ID GSE28475 [20,21]. The samples were collected from the superior frontal gyrus of the dorsolateral prefrontal cortex or the middle frontal gyrus, if the former was unavailable. For the current work, only the samples acquired by DASL-based assays were included in the analysis to keep the data consistent. Since the gene-expression levels are dependent on age, as shown by another study by the same group [22], we chose to use the median expressions for the samples aged between 2 and 14.
A flowchart of the methodology used in this work can be seen in Figure 1. tINIT [23] was chosen as the contextualizing algorithm, and Human1 [24] was selected as the reference model, since it unifies all the information from the other three generic human metabolic models, namely HMR 2.0 [25], iHsa [26] and Recon3D [27], and has more metabolites and reactions present than the other three models. The decision for threshold selection in the tINIT algorithm was made according to the distribution of expression values. The distribution of the quantile normalized gene expression values is bimodal, with a local minimum between 9 and 9.25. We selected 9.25 as the threshold since lower values tend to increase the number of reactions and metabolites in the model while decreasing the specificity of the model ( Figure S1).
Biomedicines 2023, 10, x FOR PEER REVIEW 3 of 18 A flowchart of the methodology used in this work can be seen in Figure 1. tINIT [23] was chosen as the contextualizing algorithm, and Human1 [24] was selected as the reference model, since it unifies all the information from the other three generic human metabolic models, namely HMR 2.0 [25], iHsa [26] and Recon3D [27], and has more metabolites and reactions present than the other three models. The decision for threshold selection in the tINIT algorithm was made according to the distribution of expression values. The distribution of the quantile normalized gene expression values is bimodal, with a local minimum between 9 and 9.25. We selected 9.25 as the threshold since lower values tend to increase the number of reactions and metabolites in the model while decreasing the specificity of the model ( Figure S1). The RAVEN toolbox [28] was used during GEM contextualization, and the simulations were performed using the COBRA toolbox [29]. All related operations were conducted in MATLAB R2017b, and Gurobi optimization software (version 8.1.1) [30] was used as the linear programming solver. In both the ASD-specific and neurotypical models, we checked if the pruned models contained all the reactions regarding the glutamineglutamate-GABA cycle, and if absent, we added them manually. Additionally, if the exchange and transport reactions for these and other essential metabolites are missing from the pruned models but not from the reference Human1, they were included in the models.

GEM Constraints
For the ASD and control models, separate reaction bounds were set if the information was available. The glucose uptake rate was set to the cerebral metabolic rate (CMR) of glucose in the frontal cortexes of the ASD and neurotypical children [31]. The maximum oxygen uptake rates were set according to the metabolic ratios of CMRO2/CMRglu for the young ASD patients and controls [32], and the minimum oxygen uptake rates were set to the minimum of CMRglu in the mild asphyxia infants [33]. The energy expenditure is estimated to be nearly 10 µmol ATP/(g min) for non-signaling cellular activities and 30 µmol ATP/(g min) for signaling-related processes, and these were set as lower (non-signaling) and upper (signaling + non-signaling) boundaries for the non-growth-associated ATP maintenance reaction [34]. The maximum reaction rates for the glutaminase [35], glutamine synthetase [36], and glutamate carboxylase [37] reactions were set as the greater Vmax value between astrocytes and neurons. For the exchange reactions of other essential nutrients, such as amino acids, the cerebrospinal fluid (CSF) concentrations were converted to flux units using CSF flow rates, and these were assumed to be equal to the maximum uptake rates for the corresponding reactions [38].

Model Simulations
Most steady-state analysis methods, such as flux balance analysis (FBA), require a predefined objective function for the optimization of the problem. For organisms like bacteria, the selection of the objective function is mostly straightforward; they aim to The RAVEN toolbox [28] was used during GEM contextualization, and the simulations were performed using the COBRA toolbox [29]. All related operations were conducted in MATLAB R2017b, and Gurobi optimization software (version 8.1.1) [30] was used as the linear programming solver. In both the ASD-specific and neurotypical models, we checked if the pruned models contained all the reactions regarding the glutamine-glutamate-GABA cycle, and if absent, we added them manually. Additionally, if the exchange and transport reactions for these and other essential metabolites are missing from the pruned models but not from the reference Human1, they were included in the models.

GEM Constraints
For the ASD and control models, separate reaction bounds were set if the information was available. The glucose uptake rate was set to the cerebral metabolic rate (CMR) of glucose in the frontal cortexes of the ASD and neurotypical children [31]. The maximum oxygen uptake rates were set according to the metabolic ratios of CMRO 2 /CMRglu for the young ASD patients and controls [32], and the minimum oxygen uptake rates were set to the minimum of CMRglu in the mild asphyxia infants [33]. The energy expenditure is estimated to be nearly 10 µmol ATP/(g min) for non-signaling cellular activities and 30 µmol ATP/(g min) for signaling-related processes, and these were set as lower (nonsignaling) and upper (signaling + non-signaling) boundaries for the non-growth-associated ATP maintenance reaction [34]. The maximum reaction rates for the glutaminase [35], glutamine synthetase [36], and glutamate carboxylase [37] reactions were set as the greater V max value between astrocytes and neurons. For the exchange reactions of other essential nutrients, such as amino acids, the cerebrospinal fluid (CSF) concentrations were converted to flux units using CSF flow rates, and these were assumed to be equal to the maximum uptake rates for the corresponding reactions [38].

Model Simulations
Most steady-state analysis methods, such as flux balance analysis (FBA), require a predefined objective function for the optimization of the problem. For organisms like bacteria, the selection of the objective function is mostly straightforward; they aim to maximize their growth and proliferate. However, for human cells, especially for brain cortex cells, maximization of the growth function is not applicable. Thus, we maximized the reaction rates for the glutaminase, glutamine synthetase, and glutamate carboxylase reactions to provide the replenishment of the glutamine in the neurotransmitter pool through the glutamine-glutamate-GABA cycle. Uniform sampling was used in order to obtain the distribution of attainable flux values rather than one of the optima obtained by the FBA. The lower and upper bounds of the reactions were set to the minimum and maximum flux values obtained through flux variability analysis (FVA) with a loopless solution option by maximizing the glutamine-glutamate-GABA cycle, if possible. A total of 10,000 samples were collected using the artificial centering hit-and-run (ACHR) algorithm within the COBRA Toolbox.

Statistical Analysis
All statistical analyses and data visualizations were performed in R software version 3.6.3. Due to the large dataset, effect size estimations were reported instead of the significance of the difference between group means/medians. The effect sizes between samples were determined using Cliff's d using the effsize R package [39]. The Cliff's d values were computed using an unbiased estimate formula, non-normal distribution assumption, and 95% confidence interval levels.

Model, Constraint, and Simulation Features
The selection of the flux constraints for model reactions is a significant step in genomescale modeling. The phenotype of an organism or cell is determined by its environment and thermodynamic limitations, as well as by its genotype. Thus, appropriate constraints should be set for the simulations to obtain a realistic set of solutions. Since the metabolisms of children and adults are dissimilar, the search for the constraints was concentrated on experimental data for children aged 2-14, but if the information was not available, the values obtained from the adult ASD and control cases were used.
The extracellular fluid of the brain and spinal cord consists of blood plasma, cerebrospinal fluid (CSF), and interstitial fluid (ISF). The brain tissue is surrounded by ISF, while the CSF fills the brain ventricles and subarachnoid space. The compositions of ISF and CSF are very similar, although they are significantly different from the composition of blood since the blood-brain barrier is a highly selective barrier that allows the passage of small non-polar molecules and lipophilic molecules and blocks the passage of molecules such as neurotransmitters out of the brain [40]. The compositions of these extracellular fluids are substantial for the brain's functioning. One study showed that the concentrations of several amino acids in CSF, such as glycine, aspartate, and arginine, are significantly different in ASD children compared to other patients with neurological problems [38]. The number of significantly different amino acid concentrations might be even higher if ASD and neurotypical children are compared; however, since it is not ethical to perform a lumbar puncture on healthy children to obtain CSF, no studies have investigated this difference. Yet, the presented information may give insight into the metabolism differences between ASD patients and others. Thus, the maximum amino acid uptake rates are set to their concentrations in CSF after the necessary conversions are done.
The glucose utilization rates were also determined to be not significantly different between ASD and control children aged between 2 and 18 [31]. The metabolic ratio of the cerebral metabolic rate of oxygen to that of glucose has not been reported for ASD; thus, in this study, the value for young adults was used to determine the maximum uptake rate of the oxygen.
The human prefrontal cortex mostly consists of astrocytes, glutamatergic excitatory pyramidal neurons, and GABAergic inhibitory interneurons. A simplified illustration of the exchanges between the astrocytes and neurons can be seen in Figure 2. Glutamate and gamma-aminobutyric acid (GABA) are two essential neurotransmitters: the former acts as excitatory, and the latter acts as an inhibitor. Neurons lack the ability to produce glutamate de novo from glucose since pyruvate carboxylase is not present in neurons [41]. Glutamate is converted into glutamine in astrocytes and released into extracellular space to later be used as precursors for glutamate and GABA in glutamatergic and GABAergic neurons, respectively. Since the neurotransmitter exchange between these cells is crucial for the functioning of the prefrontal cortex, the objective functions for both models are here selected as the maximization of the glutaminase, glutamine synthetase, and glutamate carboxylase reactions for the flux variability analysis, which is later used as the flux bounds of the model reactions for the uniform sampling analysis. Although the ACHR sampling algorithm used here randomly picks an objective function among the model's reactions for each sampling, constraining all reactions to the FVA results provided non-zero flux through the glutamine-glutamate-GABA cycle for each sample.
gamma-aminobutyric acid (GABA) are two essential neurotransmitters: the former acts as excitatory, and the latter acts as an inhibitor. Neurons lack the ability to produce glutamate de novo from glucose since pyruvate carboxylase is not present in neurons [41]. Glutamate is converted into glutamine in astrocytes and released into extracellular space to later be used as precursors for glutamate and GABA in glutamatergic and GABAergic neurons, respectively. Since the neurotransmitter exchange between these cells is crucial for the functioning of the prefrontal cortex, the objective functions for both models are here selected as the maximization of the glutaminase, glutamine synthetase, and glutamate carboxylase reactions for the flux variability analysis, which is later used as the flux bounds of the model reactions for the uniform sampling analysis. Although the ACHR sampling algorithm used here randomly picks an objective function among the model's reactions for each sampling, constraining all reactions to the FVA results provided nonzero flux through the glutamine-glutamate-GABA cycle for each sample.

Cellular Respiration and Energy Metabolism Is Disturbed in ASD
Mitochondria is the main ATP generation unit of the cells through oxidative phosphorylation and is also responsible for the regulation of reactive oxygen species (ROS). Although only 5% of ASD children are diagnosed with mitochondrial disease, up to 80% of ASD children show signs of mitochondrial dysfunction [16]. Mitochondrial dysfunction biomarkers were found to be elevated in the blood, plasma, and urine samples of individuals with ASD [42]. On the other hand, a large amount of cellular ROS is produced in mitochondria as a byproduct of electron transport chain (ETC) activity. The brain only constitutes 2% of a human's body weight, yet it accounts for nearly 20% of the body's required daily energy. Glucose is the primary energy source for brain cells, and glycolysis, the tricarboxylic acid (TCA) cycle and oxidative phosphorylation are the mechanisms for producing ATP by breaking glucose. Theoretically, a maximum of 36 ATPs can be produced per consumed glucose molecule, two produced by glycolysis, two by the TCA cycle, and 32 by oxidative phosphorylation. Several studies detected that the

Cellular Respiration and Energy Metabolism Is Disturbed in ASD
Mitochondria is the main ATP generation unit of the cells through oxidative phosphorylation and is also responsible for the regulation of reactive oxygen species (ROS). Although only 5% of ASD children are diagnosed with mitochondrial disease, up to 80% of ASD children show signs of mitochondrial dysfunction [16]. Mitochondrial dysfunction biomarkers were found to be elevated in the blood, plasma, and urine samples of individuals with ASD [42]. On the other hand, a large amount of cellular ROS is produced in mitochondria as a byproduct of electron transport chain (ETC) activity. The brain only constitutes 2% of a human's body weight, yet it accounts for nearly 20% of the body's required daily energy. Glucose is the primary energy source for brain cells, and glycolysis, the tricarboxylic acid (TCA) cycle and oxidative phosphorylation are the mechanisms for producing ATP by breaking glucose. Theoretically, a maximum of 36 ATPs can be produced per consumed glucose molecule, two produced by glycolysis, two by the TCA cycle, and 32 by oxidative phosphorylation. Several studies detected that the genes that play a role in oxidative phosphorylation are downregulated in the brains of autistic patients [43,44].
In our genome-scale metabolic brain model, the attainable flux distributions for the aerobic respiration-related pathways showed that the glycolysis pathway is more active in the ASD model (Table 1 and Figure 3). In parallel with these results, one study conducted on lymphoblastoid cell lines (LCLs) from ASD patients, their siblings, and controls demonstrated that the glycolytic capacities of ASD LCLs were 19% and 12% greater than their siblings and unrelated controls, respectively [45]. While the ratios of the average of these calculated fluxes are similar through glycolysis, the ratios are variable through the TCA cycle. The fluxes converting citrate, cis-aconitate, and succinate are faster, and fluxes converting isocitrate, α-ketoglutarate, malate, and oxaloacetate are predicted to be at least 40% slower in the ASD model. The inconsistent flux levels observed in the TCA cycle of the ASD model compared to the control model suggest that metabolites downstream of isocitrate are differentially produced in ASD and that there is a disruption in the TCA cycle. One may speculate that two of these consequent downstream enzymes, isocitrate dehydrogenase (IDH) and α-ketoglutarate dehydrogenase (α-KGDH), are activated and regulated by calcium ions (Ca 2+ ) in the mitochondrial matrix [46], and the deficiency of Ca 2+ may create disruptions in the mitochondrial metabolism. Indeed, genetic variants of voltage-gated calcium channels (VGCCs), transmembrane proteins mediating the Ca 2+ flux to neurons, were detected in ASD and speculated to be a part of the pathophysiology of ASD [47]. Alterations in TCA metabolite levels were also detected in other studies [48][49][50]. The changes in the TCA metabolites may indicate that the TCA cycle is a common pathway impacted in ASD. Consequently, as our computational findings point out, the increase in glycolysis rates in ASD may be related to the disruption in the TCA cycle. As less ATP is produced through the TCA cycle, an increase in the glycolytic pathway may be needed to meet the ATP demand of the cell. The anomaly in the function of the TCA cycle can also be deducted from the ratio of CO 2 secreted for every glucose molecule uptook. Theoretically, 6 CO 2 molecules are formed per glucose molecule by complete aerobic respiration, while none are formed after glycolysis. As seen in Figure 4, the CO 2 secretion per glucose uptake for the ASD-specific model is found to be lower compared to the control model with a large effect size, which is an indicator of incomplete aerobic respiration. The disturbance of the TCA cycle reduces the aerobic respiration rate, and this phenomenon generally results in elevated levels of pyruvate, as well as its derivatives, lactate and alanine [16]. In our metabolic brain model, the secretion of pyruvate is found to be significantly greater in the ASD model and lactate is found to be generally transferred into  Table 1. ASD-specific metabolic traits predicted by the genome scale metabolic brain model.

ASD-Related Changes in Mitochondrial Dysfunction and Oxidative Stress Related Markers
The mitochondrial electron transport chain (ETC) is located in the inner membrane of the mitochondria and produces ATP by generating a proton gradient between the mitochondrial matrix to the inner membrane space. ETC consists of 5 complexes: complex I (NADH dehydrogenase), complex II (succinate dehydrogenase), complex III (cytochrome bc1 complex), complex IV (cytochrome c oxidase), and complex V (ATP synthase) ( Figure 5). We computationally observed some unusual results ( Figure 5  conducted on brain tissues, and ETC complex functions can vary in different tissues and even in the same tissue. Oxidative stress is believed to be a chronic condition in ASD as no studies have found a correlation between age and the levels of oxidative stress markers in ASD patients [59]. Glutathione peroxidase (GPx) is a significant antioxidant enzyme that takes part in the peroxyl scavenging mechanism by converting GSH and hydrogen peroxide (H2O2) to GSSG and water. Our results reveal that the total rate of GPx-catalyzed reactions in the cytosol and mitochondria have similar flux ranges in both models, indicating that both the ASD and control models scavenge similar amounts of H2O2 (Figure 4). In the literature, the GPx activities in the erythrocytes and plasma of ASD patients were found to be low [60,61], as well as in the cerebellum [62], with two conflicting results [63,64] suggesting there is not a pattern regarding GPx in ASD. The reaction converting GSSG to GSH is catalyzed by glutathione-disulfide reductase (GSR), which is considered to play a significant part in antioxidant production in cells. Our simulation results show that the fluxes for the reaction catalyzed by GSR are significantly increased in ASD cases, suggesting inflated oxidative stress in the prefrontal cortex of ASD individuals. However, the role of The disturbance of the TCA cycle reduces the aerobic respiration rate, and this phenomenon generally results in elevated levels of pyruvate, as well as its derivatives, lactate and alanine [16]. In our metabolic brain model, the secretion of pyruvate is found to be significantly greater in the ASD model and lactate is found to be generally transferred into the cell, yet the ASD model needed less lactate than the control model. Pyruvate is converted to alanine by alanine aminotransferase (ALT), whereas the pyruvate-derivative oxaloacetate is converted to aspartate by the aspartate aminotransferase (AST) enzyme. Our simulations show that pyruvate secretion and the activity of reactions catalyzed by ALT and AST are higher in the ASD model, thus supporting the reduced aerobic respiration hypothesis ( Figure 4). Several studies also revealed elevated ALT and AST levels in the bodily fluids of ASD patients [51], confirming our results. The disturbances in the TCA cycle may also be explained by the elevated levels of ammonia in the cell, which is a byproduct of protein metabolism and is accepted as another biomarker of the mitochondrial disorder [52][53][54]. Our results showed that the ASD model secretes higher ammonia on average (Figure 4). Elevated ammonia levels may constrain the TCA cycle in neurons and glia by inhibiting cycle enzymes such as pyruvate dehydrogenase (PDH), isocitrate dehydrogenase, (IDH), and α-ketoglutarate dehydrogenase (α-KGDH).

ASD-Related Changes in Mitochondrial Dysfunction and Oxidative Stress Related Markers
The mitochondrial electron transport chain (ETC) is located in the inner membrane of the mitochondria and produces ATP by generating a proton gradient between the mitochondrial matrix to the inner membrane space. ETC consists of 5 complexes: complex I (NADH dehydrogenase), complex II (succinate dehydrogenase), complex III (cytochrome bc1 complex), complex IV (cytochrome c oxidase), and complex V (ATP synthase) ( Figure 5). We computationally observed some unusual results ( Figure 5). Among five complexes, ETC complexes I and V show higher flux values in the ASD model, while the remaining complexes have lower fluxes, with all complexes showing large Cliff's d values. The ETC complex protein levels were measured in the cerebellum, frontal, parietal, occipital, and temporal cortices of ASD patients and controls of different ages, and these ETC complex levels were found to be lower in the cerebellum, frontal cortex, and temporal cortex of the ASD cases aged between 4 and 10 compared to the age-matched controls, but not in other cortices and ASD cases aged between 14 and39 [55]. Specifically, the levels of complex I are significantly decreased in the frontal cortex of autistic children, with non-significant decreases in levels of other complexes. Another study by the same group compared the ETC activities in the frontal cortex of ASD patients and detected that 9 out of 14 ASD patients had at least one ETC complex with lower activity than controls. The most affected complexes were found as complexes I and V [56]. The decrease in the complex I activity in blood cells was also been found in recent studies [51]. The activities for all the ETC complexes in the ASD granulocytes were also significantly lower [57]. Another late study also detected abnormal complex activities in the ASD fibroblasts [58]. Although our simulation results on the ETC complex ( Figure 5) were unanticipated, none of the aforementioned activity studies were conducted on brain tissues, and ETC complex functions can vary in different tissues and even in the same tissue.
Oxidative stress is believed to be a chronic condition in ASD as no studies have found a correlation between age and the levels of oxidative stress markers in ASD patients [59]. Glutathione peroxidase (GPx) is a significant antioxidant enzyme that takes part in the peroxyl scavenging mechanism by converting GSH and hydrogen peroxide (H 2 O 2 ) to GSSG and water. Our results reveal that the total rate of GPx-catalyzed reactions in the cytosol and mitochondria have similar flux ranges in both models, indicating that both the ASD and control models scavenge similar amounts of H 2 O 2 ( Figure 4). In the literature, the GPx activities in the erythrocytes and plasma of ASD patients were found to be low [60,61], as well as in the cerebellum [62], with two conflicting results [63,64] suggesting there is not a pattern regarding GPx in ASD. The reaction converting GSSG to GSH is catalyzed by glutathione-disulfide reductase (GSR), which is considered to play a significant part in antioxidant production in cells. Our simulation results show that the fluxes for the reaction catalyzed by GSR are significantly increased in ASD cases, suggesting inflated oxidative stress in the prefrontal cortex of ASD individuals. However, the role of GSR in ASD is still poorly understood, with previous studies finding conflicting results. One study detected no significant difference between the ASD and control groups, yet 60% of the ASD cases showed lower and higher GSR activities than 95% CI (confidence interval) of the control group [62]. The other study detected higher GSR activity in the plasma of ASD patients [64].
Superoxide dismutase (SOD) and catalase (CAT) are two other significant antioxidant enzymes. SOD catalyzes the dismutation of superoxide anion (O − 2 ) to H 2 O 2 , whereas CAT catalyzes the conversion of H 2 O 2 to H 2 O. Our simulation result for the SOD enzyme is displayed in Figure 4. The mean of the reaction flux is lower in the ASD model, but the range for the flux is wider than that in the control model, which suggests that the ASD model has a higher variability for SOD activity compared to the control model in accordance with several literature findings. In fact, literature studies have presented conflicting results. Lower [61], higher [65][66][67], or similar [63,68] SOD activities in ASD plasma and erythrocytes have been reported in comparison to controls. The same trend was also detected for CAT: lower [66,67], higher [65,69], and non-significant difference [68]. The mean values of CAT activities are closer in our ASD and control models, with a slight elevation in ASD, yet the maximum CAT activities are greater in the ASD model. The overactivity of oxidative stress-protective enzymes such as CAT, GSR, and SOD in our ASD model may indicate that these enzymes should overwork to compensate for the oxidative stress in the ASD brain.

Neuroinflammatory Markers in ASD
Neuroinflammation is the inflammation of the central nervous system and is indicated to take part in the pathogenesis of several neurological diseases. The activation of glial cells, such as astrocytes, which are one of the most abundant cells in the prefrontal cortex, is associated with the neuroinflammation status of the brain and releases pro-inflammatory mediators. The inflammatory status of ASD children has been indicated by the increase in the pro-inflammatory cytokines [12]. In this section, we present our computational findings on exchange fluxes of neuroinflammation-related compounds and also elaborate on the metabolite biomarkers for inflammation.
Inositol (also known as myo-inositol) and choline-containing compounds are two indicators of neuroinflammation in the brain and are generally detected via proton magnetic resonance spectroscopy analyses. Activation of glial cells is thought to be correlated with higher inositol and choline-containing compound levels [70]. The activation in these cells is also associated with the neuroinflammation status of the brain and the release of proinflammatory mediators [71], suggesting that inositol and choline levels may be relevant to the neuroinflammation status. Our simulation results are shown in Figure 6. Acetylcholine and choline secretions are found to be higher in the ASD model with large effect sizes. Although the difference is not significant for inositol, the control model secretes more inositol on average. We may speculate that although the literature suggests no significant difference for these metabolites in prefrontal cortex [72], the neuroinflammation in the prefrontal cortex stimulates the astrocyte activation and its relevant reactions, which in turn produces more choline-related metabolites.
Phospholipases, enzymes hydrolyzing phospholipids to fatty acids, are at the center of the fatty acid mechanism. Specifically, phospholipase A2 (PLA2) breaks the phospholipids into arachidonic acid and polyunsaturated free fatty acids (PUFA), which are further hydrolyzed to eicosanoids such as prostaglandins and leukotrienes. We found that the reaction catalyzed by the PLA2 enzyme has significantly higher fluxes in the control model ( Figure 6). A small number of studies investigated the activity of PLA2 in ASD and found conflicting results; in particular, one study detected lower PLA2 activity in the serum of ASD children [73], and two others concluded higher activities in the RBC [74] and blood [75] of ASD children. Docosahexaenoic acid (DHA) is the most abundant n-3 PUFA in the brain and regulates the function of glial cells and neurons. It also protects the neurons from damage in several brain diseases and diminishes neuroinflammation [76,77]. The precursor of DHA, eicosapentaenoic acid (EPA), is also effective against inflammation, although it is less abundant than DHA in the brain [78]. Our simulations show that both DHA and EPA secretions are higher in the control model, although the difference is not very explicit, indicating that ASD people may be more prone to inflammation and cell damage ( Figure 6). Studies have shown that the DHA levels are lower in the blood of ASD children, whereas only one study reported significantly lower levels of EPA in the ASD brain [78]. The differences in these fatty acids' secretion capacities indeed indicate perturbations in fatty acid metabolism. Phospholipases, enzymes hydrolyzing phospholipids to fatty acids, are at the center of the fatty acid mechanism. Specifically, phospholipase A2 (PLA2) breaks the phospholipids into arachidonic acid and polyunsaturated free fatty acids (PUFA), which are further hydrolyzed to eicosanoids such as prostaglandins and leukotrienes. We found that the reaction catalyzed by the PLA2 enzyme has significantly higher fluxes in the control model ( Figure 6). A small number of studies investigated the activity of PLA2 in ASD and found conflicting results; in particular, one study detected lower PLA2 activity in the serum of ASD children [73], and two others concluded higher activities in the RBC [74] and blood [75] of ASD children. Docosahexaenoic acid (DHA) is the most abundant n-3 PUFA in the brain and regulates the function of glial cells and neurons. It also protects the neurons from damage in several brain diseases and diminishes neuroinflammation [76,77]. The precursor of DHA, eicosapentaenoic acid (EPA), is also effective against inflammation, although it is less abundant than DHA in the brain [78]. Our simulations show that both DHA and EPA secretions are higher in the control model, although the difference is not very explicit, indicating that ASD people may be more prone to inflammation and cell damage ( Figure 6). Studies have shown that the DHA levels are lower in the blood of ASD children, whereas only one study reported significantly lower levels of EPA in the ASD brain [78]. The differences in these fatty acids' secretion capacities indeed indicate perturbations in fatty acid metabolism.
Prostaglandins and leukotrienes are two important classes of eicosanoids present in inflammatory conditions. Prostaglandin E2 (PGE2) is one of the most abundant prostaglandins in the human body and is a mediator for many biological functions, including inflammation [79]. PGE2 may induce both pro-inflammatory and anti-inflammatory processes depending on the context, but mostly in the pro-inflammatory direction [80]. On the other hand, prostaglandin E1 (PGE1) is thought to function as an anti-inflammatory mediator for several animals, including humans [81,82]. According to our simulations, the Prostaglandins and leukotrienes are two important classes of eicosanoids present in inflammatory conditions. Prostaglandin E2 (PGE2) is one of the most abundant prostaglandins in the human body and is a mediator for many biological functions, including inflammation [79]. PGE2 may induce both pro-inflammatory and anti-inflammatory processes depending on the context, but mostly in the pro-inflammatory direction [80]. On the other hand, prostaglandin E1 (PGE1) is thought to function as an anti-inflammatory mediator for several animals, including humans [81,82]. According to our simulations, the secretion of PGE2 from the cell, often a pro-inflammatory metabolite, is higher in the ASD model, whereas the PGE1 secretion is lower. Leukotrienes, especially leukotriene B4, are indicated as pro-inflammatory mediators and suggested as potential therapeutic targets for the modulation of inflammation [83]. The leukotriene levels were found to be higher in the ASD model compared to those in neurotypical controls ( Figure 6). These results, combined with the literature findings [84][85][86][87] on the bodily fluids of ASD patients, imply that there is inflammatory stress in the prefrontal cortex of ASD individuals.

Sphingolipid Metabolism Changes in ASD
Compared to other tissues, the nervous system, especially the brain, has the highest lipid content and complexity. Sphingolipids, being one of the major classes of lipids that are essential for eukaryotic cells, are highly enriched in the brain. They have several structural roles in the plasma membranes and are also signaling molecules regulating cellular events such as cell growth, differentiation, senescence, and apoptosis [88,89]. Sphingolipids are implicated as a critical component for brain development and function. Altered sphingolipid metabolism has been implied in several neurological and psychiatric disorders. Thus far, no detailed studies have focused on the altered metabolism of sphingolipids in the brains of ASD patients and its effects on the pathophysiology of the disease. Currently, only two sphingolipid metabolites have been reported in the literature for ASD. Sphingosine-1-phosphate (S1P) was found to be significantly increased in ASD serum [90], and several ceramide types have been elevated in the ASD brain [91] and prefrontal cortex [92].
Ceramides participate in important signaling processes, such as inflammation, and can be generated by pro-inflammatory cytokines [93]. Ceramide can be phosphorylated to ceramide-1-phosphate, which may activate PLA2, and may thus induce inflammation. Ceramides can also be degraded to sphingosines, which can be converted to S1P. S1P has many roles in the cell, including cell migration, proliferation, cellular architecture, and inflammation [94], and is known to have anti-inflammatory effects [95]. Our simulations have detected several differences in the sphingolipid metabolisms of ASD and neurotypical people (Figure 7). On average, the sphingolipid metabolism has higher fluxes in the ASD model, but the ratios through the pathway are consistent with each other, except for ceramide conversion to sphingosine. The fluxes through the pathway are nearly 1.2 to 1.4 times faster in ASD, except for the conversion of ceramide to sphingosine. The S1P secretion fluxes are increased in our simulations (Figure 7), in parallel with the literature results [90], while the pro-inflammatory ceramide secretion is restrained in ASD individuals on average, contrary to our expectations. Although more ceramide is produced in the ASD model, most of the produced ceramide is converted to other derivatives like sphingomyelin or sphingosine, which may explain the decrease in ceramide secretion. The increase in ceramide production can be induced by oxidative stress [96], which is shown to be elevated in ASD children by this study. Another derivative of ceramides, glucosylceramides, are indicated as both pro-and anti-inflammatory mediators and are also found to be higher in ASD patients.

Conclusions
Our simulations on reconstructed ASD-specific and neurotypical control GEMs demonstrated that several pathways, which have already been shown to be disturbed in ASD, generally yield complementary results with literature findings. The energy metabo- On a special note, our simulations indicated altered sphingolipid metabolism between the ASD and neurotypical control models and this metabolism needs further research to shed light on the effect of sphingolipids on the pathogenesis of ASD. The research on the relationship between sphingolipid metabolism and ASD has been limited. In fact, the majority of studies on ASD omics have focused on transcriptomics or metabolomics, and it is only recently that lipidomics analyses have been conducted, although lipids constitute 36-66% of the dry weight of the human brain depending on grey or white matter, and nearly one-third of the total lipids are sphingolipids [97]. The increase in the number of lipidomics studies may enhance our understanding of the brain structure and function in ASD children.

Conclusions
Our simulations on reconstructed ASD-specific and neurotypical control GEMs demonstrated that several pathways, which have already been shown to be disturbed in ASD, generally yield complementary results with literature findings. The energy metabolism is altered in ASD, as indicated by the distinct patterns in glycolysis and the TCA cycle. One of the most prevalent comorbidities of ASD, mitochondrial dysfunction, is also supported by the disparities in the activities of ETC complexes and the increase in the mitochondrial dysfunction biomarkers. In addition, the reactions catalyzed by the oxidative stress marker enzymes, as well as the ones catalyzed by oxidative-stress protective enzymes, showed higher fluxes in ASD, indicating the presence of oxidative stress in ASD. Pro-inflammatory markers are found to have higher fluxes in ASD, whereas contradictory results are found for anti-inflammatory metabolites. Sphingolipid metabolites, especially S1P and ceramide, are promising biomarkers for ASD and should be thoroughly investigated by in vitro experiments.
Omics data, such as transcriptomics, proteomics, and metabolomics, enhance the predictive power of systems biology tools such as genome-scale modeling. However, the lack of adequate knowledge about ASD was the most significant challenge encountered throughout this modeling study. Despite the limited number of omics studies conducted on the brains of ASD patients, all of them are post-mortem studies due to inaccessibility to the area, which may have distinctive features compared to ante-mortem analyses. Additionally, the number of studies reporting the enzyme activities in the brain, specifically the prefrontal cortex, is quite low, and in most cases, they conflict with each other. Thus, the neuromodeling studies, including ours, have inadequate resources to validate the computational results. Furthermore, the heterogeneity of the disease and the complexity of the nervous system hinders the efforts to understand the underlying conditions behind ASD and the impact of altered pathways on the symptoms of ASD. However, customized forecasts of patients' phenotypic characterization can be improved by employing genome-scale metabolic models. The metabolic models tailored to individual conditions, tissues, and patients shed light on cell-specific changes in the metabolic network, allowing for the prediction of disease metabolism biomarkers. All in all, we may conclude that genomescale modeling may be beneficial in investigating the disturbed mechanisms in ASD.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biomedicines11020583/s1, Figure S1: The distribution of the median quantile normalized gene expression values for samples aged between 2 and 14. Data Availability Statement: Publicly available datasets were analyzed in this study. The constructed models for this study are available at https://doi.org/10.5281/zenodo.7602974.